﻿using System;
using System.Collections.Generic;
using System.Drawing;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace XChip_Common
{
    /// <summary>
    /// Bezier Spline methods
    /// </summary>
    public static class BezierSpline
    {
        /// <summary>
        /// 获取开放式Bezier样条曲线控制点。
        /// </summary>
        /// <param name="knots">输入结贝塞尔样条曲线点。</param>
        /// <param name="firstControlPoints">输出第一个控制点数组结。长度-1长度。</param>
        /// <param name="secondControlPoints">输出第二个控制点数组结。长度-1长度。</param>
        /// <exception cref="ArgumentNullException"><paramref name="knots"/>
        ///  参数不能为空。</exception>
        /// <exception cref="ArgumentException"><paramref name="knots"/>
        /// 数组必须至少包含两个点。</exception>
        public static void GetCurveControlPoints(Point[] knots,
            out Point[] firstControlPoints, out Point[] secondControlPoints)
        {
            if (knots == null)
                throw new ArgumentNullException("knots");
            int n = knots.Length - 1;
            if (n < 1)
                throw new ArgumentException
                ("At least two knot points required", "knots");
            if (n == 1)
            { // Special case: Bezier curve should be a straight line.
                firstControlPoints = new Point[1];
                // 3P1 = 2P0 + P3
                firstControlPoints[0].X = (2 * knots[0].X + knots[1].X) / 3;
                firstControlPoints[0].Y = (2 * knots[0].Y + knots[1].Y) / 3;

                secondControlPoints = new Point[1];
                // P2 = 2P1 – P0
                secondControlPoints[0].X = 2 *
                    firstControlPoints[0].X - knots[0].X;
                secondControlPoints[0].Y = 2 *
                    firstControlPoints[0].Y - knots[0].Y;
                return;
            }

            // Calculate first Bezier control points
            // Right hand side vector
            double[] rhs = new double[n];

            // Set right hand side X values
            for (int i = 1; i < n - 1; ++i)
                rhs[i] = 4 * knots[i].X + 2 * knots[i + 1].X;
            rhs[0] = knots[0].X + 2 * knots[1].X;
            rhs[n - 1] = (8 * knots[n - 1].X + knots[n].X) / 2.0;
            // Get first control points X-values
            double[] x = GetFirstControlPoints(rhs);

            // Set right hand side Y values
            for (int i = 1; i < n - 1; ++i)
                rhs[i] = 4 * knots[i].Y + 2 * knots[i + 1].Y;
            rhs[0] = knots[0].Y + 2 * knots[1].Y;
            rhs[n - 1] = (8 * knots[n - 1].Y + knots[n].Y) / 2.0;
            // Get first control points Y-values
            double[] y = GetFirstControlPoints(rhs);

            // Fill output arrays.
            firstControlPoints = new Point[n];
            secondControlPoints = new Point[n];
            for (int i = 0; i < n; ++i)
            {
                // First control point
                firstControlPoints[i] = new Point((int)x[i], (int)y[i]);
                // Second control point
                if (i < n - 1)
                    secondControlPoints[i] = new Point((int)(2 * knots[i + 1].X - x[i + 1]), (int)(2 * knots[i + 1].Y - y[i + 1]));
                else
                    secondControlPoints[i] = new Point(((int)(knots[n].X + x[n - 1]) / 2), (int)((knots[n].Y + y[n - 1]) / 2));
            }
        }

        /// <summary>
        /// 为第一个Bezier控制点的坐标（x或y）求解三对角系统。
        /// </summary>
        /// <param name="rhs">Right hand side vector.</param>
        /// <returns>Solution vector.</returns>
        private static double[] GetFirstControlPoints(double[] rhs)
        {
            int n = rhs.Length;
            double[] x = new double[n]; // Solution vector.
            double[] tmp = new double[n]; // Temp workspace.

            double b = 2.0;
            x[0] = rhs[0] / b;
            for (int i = 1; i < n; i++) // Decomposition and forward substitution.
            {
                tmp[i] = 1 / b;
                b = (i < n - 1 ? 4.0 : 3.5) - tmp[i];
                x[i] = (rhs[i] - x[i - 1]) / b;
            }
            for (int i = 1; i < n; i++)
                x[n - i - 1] -= tmp[n - i] * x[n - i]; // Backsubstitution.

            return x;
        }
    }
}
